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The inability of otherwise successful dynamical models to reproduce the "HBT radii" extracted 
from two-particle correlations measured at the Relativistic Heavy Ion Collider (RHIC) is known as 
the "RHIC HBT Puzzle". Most comparisons between models and experiment exploit the fact that 
for Gaussian sources the HBT radii agree with certain combinations of the space-time widths of 
the source which can be directly computed from the emission function, without having to evaluate, 
at significant expense, the two-particle correlation function. We here study the validity of this 
approach for realistic emission function models some of which exhibit significant deviations from 
simple Gaussian behaviour. By Fourier transforming the emission function we compute the 2-particle 
correlation function and fit it with a Gaussian to partially mimic the procedure used for measured 
correlation functions. We describe a novel algorithm to perform this Gaussian fit analytically. 
We find that for realistic hydrodynamic models the HBT radii extracted from this procedure agree 
better with the data than the values previously extracted from the space-time widths of the emission 
function. Although serious discrepancies between the calculated and measured HBT radii remain, 
we show that a more "apples-to-apples" comparison of models with data can play an important role 
in any eventually successful theoretical description of RHIC HBT data. 

PACS numbers: 25.75.Ld, 25.75.Gz, 24.10.Nz, 25.75.-q 



I. INTRODUCTION 

Two-particle intensity interferometry is widely used to 
characterize the space-time aspects of the freeze-out con- 
figuration in relativistic heavy ion collisions 1]. It is 
common to condense this information in terms of char- 
acteristic length scales of the "homogeneity regions" 2j 
from which particles of a given momentum originate. 

In this paper we discuss the degree to which homo- 
geneity lengths extracted in quite different ways may be 
validly compared. Throughout our study, we restrict 
ourselves to interference effects between identical, non- 
interacting bosons, resulting from Bose-Einstein statis- 
tics. Since final state interactions (e.g. Coulomb effects) 
affect most interferometry studies, our study may be re- 
garded (1) as a proof-of-principle example that care must 
be taken to perform "apples-to-apples" comparisons, and 
(2) as an estimate of the magnitude of the differences for 
two popular theoretical models. 

The homogeneity length scales are extracted in exper- 
iments by assuming that the homogeneity region can be 
approximated by a Gaussian-profile ellipsoid in configu- 
ration space, resulting in a Gaussian two-particle momen- 
tum correlation function, and performing a semi-analytic 
Gaussian fit to the relative momentum dependence of the 
measured correlation function (see e.g. Il| for details). 
Following common practice, we will refer in the following 
to the size parameters obtained from Gaussian fits to the 
correlation function as "HBT radii" . 

Fitting experimental data to functional forms other 
than Gaussian is common in studies of elementary parti- 
cle collisions, for which Gaussian fits clearly fail. In heavy 
ion collisions, the Gaussian ansatz works relatively well, 
but, especially with the high quality and high-statistics 
data sets now available at RHIC, finer, non-Gaussian 



structures may be physically interesting. Instead of in- 
venting ad-hoc functional forms with which to fit the 
correlation functions, or functionally expanding about a 
Gaussian fitting form 0, Q , imaging 0, IE the homo- 
geneity region is perhaps the most promising route to 
explore these structures. In this paper we do not take 
up this issue. Instead, we note that most experimental 
studies in heavy-ion physics to date have used the Gaus- 
sian ansatz (jj, and we explore some ways in which HBT 
radii obtained in this way from data may be compared 
to model calculations. 

If the homogeneity region is indeed Gaussian in profile, 
then the HBT radii agree exactly with appropriate com- 
binations of the root-mean-squared (RMS) variances of 
its spatial distribution Q . Given a theoretical model for 
the freeze-out configuration, calculating these space-time 
variances is much easier than computing and fitting the 
correlation function. Many comparisons between mod- 
els and data therefore use this short-cut, comparing the 
space-time variances directly to the experimental HBT 
radii. However, since the homogeneity region is seldom 
perfectly Gaussian, such direct comparisons are question- 
able. 

This raises the question to what extent some of the 
persistently observed discrepancies between model pre- 
dictions and measurements of the HBT radii (the so- 
caUed "RHIC HBT Puzzle") might be due to such an 
"apples-with-oranges" comparison. Indeed, HBT radii 
calculated with Boltzmann/cascade models which are 
based on Gaussian fits to the simulated correlation func- 
tions agree somewhat better with measurements than 
do radii based on an extraction of space-time variances 
from hydrodynamic calculations Q]. Whether this is 
due to a more realistic modeling of the collision in 
the Boltzmann/cascade approach or the shortcomings 
of the comparison of variances with HBT radii in the 
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hydrodynamic case is unclear. Similarly, differences 
between hydrodynamic calculations of space-time vari- 
ances 0, and Gaussian HBT radii fitted to three- 
dimensional [TH IT^ ITsj and one-dimensional 13 Q| 
correlations have been observed. However, since these 
calculations were performed using different initial condi- 
tions and other parameters, it is unclear whether this, or 
the different extraction methods, were responsible for the 
observed differences. Here, we focus on the different ex- 
traction techniques using the same hydrodynamic model 
and parameters. 

One cascade model (MPC which reports RMS 

variances shows discrepancies with data similar to the 
hydrodynamic models. Studies 0, ITsL IT9l | performed 
within the Boltzmann/cascade framework show that 
space-time variances of the freeze-out configuration and 
Gaussian fits to the correlator can yield quite different 
radius parameters, mostly due to long tails in the spa- 
tial freeze-out distribution from resonance decays which 
strongly affect the space-time variances but are not re- 
flected by Gaussian fits to the correlation function, ac- 
cording to hydrodynamic calculations (See, how- 
ever, the recent study by Kisiel et al jU, which addresses 
this issue in detail in the context of a blast-wave param- 
eterization.) Hydrodynamic calculations of the space- 
time variances therefore usually do not include resonance 
decay contributions in the emission function Q. Still, 
the comparison in involves two differently determined 
quantities, and in the present paper we eliminate this 
shortcoming. 

To do so requires two additional steps beyond the cal- 
culation of the model emission function: (i) The correla- 
tion function must be computed via Fourier transforma- 
tion (for noninteracting identical particles) or by folding 
with a relative wave function that includes final state in- 
teraction effects (for with long-range final state interac- 
tions). This is straightforward albeit numerically expen- 
sive since it involves multiple space-time integrals, (ii) 
A Gaussian fit to the three-dimensional correlation func- 
tion must be performed, including a correlation strength 
parameter A as in the experiment. 

We here concentrate on non-interacting pairs of iden- 
tical particles as the practically most important case and 
also in order to simplify as much as possible the com- 
putation of the correlator. For the second step we de- 
velop an analytical Gaussian fit algorithm which reduces 
the multi-dimensional fit problem to a simple set of lin- 
ear equations for diagonalizing a four-dimensional ma- 
trix. This should help theoretical modelers to overcome 
the barrier of unfamiliarity when faced with a multi- 
parameter fitting problem. 

We apply our procedure to emission functions from hy- 
drodynamic calculations 9J and from the blast-wave pa- 
rameterization j2^. Both generate non-Gaussian freeze- 
out distributions, due in large measure to finite-size ef- 
fects coupled with strong collective fiow which is known 
to be important at RHIC. On the way, we also discuss 
and analyze Gaussian fits to 1-dimensional projections of 



the 3-dimensional correlator. This allows for comparison 
with earlier work along these lines [T^ |20| and first in- 
troduces our new analytic Gaussian fit algorithm in an 
easy and transparent simpler setting. 



II. VARIANCES VERSUS HBT RADII 

Experimentally, the correlation function between two 
identical particles, as a function of their relative mo- 
mentum q = Pa~Pb and their average (pair) momentum 
K={Pa+Pb)/2, is given by 



C{q,K) 



A{q,K) 
B{q,K)- 



(1) 



where A{q, K) is the signal distribution and B{q, K) is 
the reference or background distribution which is ideally 
similar to A in all respects except for the presence of fcm- 
toscopic correlations (see e.g. [l| for details). C{q,K) is 
the modification to the conditional probability for mea- 
suring particle h with momentum pi, = K—^q if particle 
a has been measured with momentum Pa = K+^q, due 
to two-particle effects sensitive to space-time separation, 
to. The explicit Jf-dependence reflects the fact that the 
separation distribution may depend on the average mo- 
mentum of the pair 2j and in general does so for explod- 
ing sources |22] . 

Theoretically, the correlation function can be calcu- 
lated from the emission function S{p, x) describing the 
probability to emit a particle from spacetime point x with 
momentum p, by convoluting it with the two-particle rel- 
ative wave function p] . For pairs of non- interacting iden- 
tical particles one has simply 0,0 



C(q,K)«l 



/ d'^xS{K,x) e"i-- 



^SxS{K,x) 



(2) 



Here q ■ x = q^t q ■ x, with = Ea—Ei, — f3 ■ q where 
/3 = K/K° = 2K/{Ea+Eb) is the average velocity of the 
pair. The « sign in Eq. (0) indicates the "smoothness ap- 
proximation" which replaces both Pa and pb by K inside 
the emission functions in the denominator 3]. Equation 
(O can be decomposed as 

C(q, K)^l + (cos(g • x)f + {sm{q ■ x))^ (3) 

where (...) indicates the (Jf-dependent) space-time av- 
erage with the emission function: 



_ Jd^xfix)SiK,x) 
~ Jd'ixS{K,x) 



(4) 



If S{K,x) is a four-dimensional Gaussian distribution 
of freeze-out points, the correlation function will likewise 
be Gaussian in the relative momentum q. It takes a par- 
ticularly simple form for midrapidity pairs (with vanish- 
ing longitudinal pair momentum, = 0) from central 
collisions between equal- mass spherical nuclei [3, @ : 



C{q)^l + Xe-(^oRl+iiRl+^fR^)_ 



(5) 



3 



Here Qo, Qs, qi are the relative momentum components 
in the Bertsch-Pratt ( "out-side-long" ) coordinate system 
0, 0|. The pair momentum dependence of the corre- 
lation function C(q, K) leads to K'-dependencies of the 
"HBT radii" i?o, Rs, and Ri (which characterize the rela- 
tive momentum widths of the correlation function) and of 
the "correlation strength" A. For fully chaotic theoretical 
Gaussian sources A = 1, but for experimental correlation 
functions usually A < 1 . Even though we here perform 
a theoretical model analysis, we keep A as a parameter 
because Gaussian fits to non-Gaussian correlation func- 
tions generally also yield A^^ 1, and experimentally such 
non-Gaussian effects on the extracted A cannot be sepa- 
rated from other origins of reduced correlation strength 
(such as contamination from misidentified particles and 
contributions from resonance decays jlj]). The HBT radii 
defined by Eq. (O convey all available geometric infor- 
mation about the source S{K,x). 

For Gaussian sources the radius parameters Rq, Rs, 
and Rl can be calculated directly from the source distri- 
bution S as RMS variances. For midrapidity pairs with 
Kl = one finds 

Rl = {xl)-2p{iJ)+^^{P), 

Rl = {il), Rf - (if), (6) 

where /3 — Kt / is the magnitude of the (transverse) 
pair velocity (which points in the Xo direction), and 

= - (a;^) (7) 

denotes the distance from the (/^-dependent) center of 
the homogeneity region for particles with momentum K. 

Experimentalists commonly extract HBT radii by fit- 
ting their experimental correlation functions ^ with the 
functional form jSjl- In contrast, most (but not all) the- 
oretical model predictions for HBT radii are based on 
a calculation of the space-time variances of the emis- 
sion function and assuming the validity of Eqs. 10 which 
holds for Gaussian sources. Of course, there is no a priori 
reason to expect a source with a perfectly Gaussian pro- 
file. Even the simplest fiow-dominated freeze-out param- 
eterizations produce clear non-Gaussian tails and edges 
[2^ . On the experimental side, high-statistics measure- 
ments show non-Gaussian behaviour, which is, however 
rarely treated quantitatively Q . In the presence of such 
non-Gaussian features, the issues are (1) whether the two 
approaches yield significantly different results, and (2) 
whether either method characterizes the physically inter- 
esting length scales of the source sufficiently well. Here, 
we address the first issue in the context of blast- wave and 
hydrodynamic models. 

Our calculations do not include experimental "noise" , 
particle mis-identification, or contributions from the de- 
cay of long-lived resonances which can reduce the fit pa- 
rameter A in Eq. ((SJ from its theoretical value of unity 
0, |2^. Instead, this parameter absorbs (and reflects) 
some of the effects of fitting a non-Gaussian function to 
a Gaussian form. This will, of course, also happen in 



experiment whenever the correlation function deviates 
from a simple Gaussian. This particular contribution to 
the fitted correlation strength A has so far received little 
attention. The model results presented here should help 
to assess the possible infiuence of non-Gaussian features 
in the data on the fitted values of A. 



III. DIRECT CALCULATION OF HBT RADII 

As explained in the Introduction, we here use model 
emission functions to compute the correlation function 
according to Eqs. H2I3II and then fit the latter with a 
Gaussian, using a procedure very similar to the one used 
in experiment. The main difference is that the theoreti- 
cal correlation function can be calculated with arbitrary 
precision, so the notion of a statistical error does not en- 
ter. Still, we will see that the fitting problem can be 
formulated in a quite analogous way. 

In the following subsection we introduce the algorithm 
for Gaussian fits through 1-dimensional cuts or projec- 
tions of the 3-dimensional correlation function. The full 
algorithm for 3-dimensional Gaussian fits is presented in 
Sec. lmBl 



A. One-dimensional Gaussian fits 

In Section VI of their paper, Wiedemann and Heinz 
|20l | calculated correlators for various model emission 
functions and extracted parameters from fits to one- 
dimensional slices of the three-dimensional correlation 
function. Although those authors called them "HBT 
radii" , we will call them "ID radii" to distinguish them 
from radii extracted from full three-dimensional fits of 
the type performed by experimentalists. 

In a given direction i (i = 0, s, I) they calculate the cor- 
relator along one of the axes i: C{qi \ qj^i=Q). They then 
find the ID radius Rfj^ i and the "directional lambda 
parameter" A^ which best approximates the correlator 
according to 

C(g,;(Z,^,=0)«l + A,e-«'^i-.'. (8) 

In particular, they calculated the correlator for a set of N 
values Qj-'^'' (similar to experimentally binning the correla- 
tion function into iV q-bins) and minimized numerically 
the quantity 

N 2 

E ^ '?5a=o) - 1) - in^'^ + Ri^A^il'^r] ■ 

k=l 

(9) 

This is reminiscent of the quantity typically minimized by 
experimenters, although in this case one also takes into 
account the experimental uncertainty of the measured 
correlator by weighting each term in the sum (bin) with 
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the inverse experimental error: 



XlD,i 
N 

E 

k=l 



(10) 



ln[C(gr;gjg.=0)-l] - In A, 
a' ('^^ 



Here, cr |p ^ represents the uncertainty in bin k on the 

quantity to be fitted, namely ln[C(gj-'^^; qj^'-^O) — l] • It 

(k) 

is related to the uncertainty aj^j-j ^ on the measured cor- 
relator C(gf ^ ^j^^^^O) itself by 



/ (fc) 

' lD,i 



d\n[C{q^ 



=0) - 1] 



'lD,i 



^z- (fc) (fc) 



=0) 



(11) 



Minimization of the quantity as in ^20] is equiva- 
lent to setting all uncertainties <j'[^ ^ to the same con- 
stant value, independent of k. However, uncertainties 
on experimental correlation functions typically have ap- 
proximately constant (fc-independent) uncertainties on 

the bin contents C(gf^gj^\=0) themselves |2|. Al- 
though statistical uncertainties on calculated correlators 
may in principle be vanishingly small, the weighting fac- 
tor [C'(g|'°''; gj^''j=0) — l] ^ which appears in Eq. H10|l as a 
result of Eq. will in general affect the resulting fit 
parameters. We choose to mimic the experimental sit- 
uation by minimizing Eq. H1Q(I . assuming constant (i.e. 
fc-independent) and infinitcsimally small errors on C, 

Minimizing Xiu i in Eq. pU|) with respect to the fit 
parameters In and i?^^ ^ by setting 



lD,i 



= 0, 



91n Ai 

we find after minimal algebra 



= 0, 



In A,; 
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where the quantities 



N 

Xn,i — 

k=\ 



-1 CT' 



Y — 



N 

E 



(fc) ^2 



(in 



-1 CT' 



(fc) ^ 

IDA) 



(12) 

(13) 
(14) 

0)-l], (15) 
(16) 



out from the ratios in Eqs. (|13I14|I . so the limit ctid,; — > 
mentioned above is well-defined. 

Minimization of Xib i differs significantly from the ex- 
perimentalists' three-dimensional fits. In particular, it 
assumes complete factorization of the correlation func- 
tion in the o, s, I directions. For at least two reasons, this 
need not be so in reality: 

(i) In a full three-dimensional fit, the three directions 
are coupled by requiring a single A parameter, inde- 
pendent of direction i. After all, according to Eq. (jHJ 
lim|q|^o C{q) = lim^^^o C {qu qj^i^o) = 1 -|- Ai should be 
independent of direction i. Thus, allowing "directional 
lambda parameters" may cause the ID fits to differ sig- 
nificantly from 3D fits. 

(ii) Perhaps more importantly, fitting separately along 
the qi axes accounts for only a set of zero measure of the 
full three-dimensional correlation function. In particu- 
lar, the correlation function may contain in the exponent 
terms such as q^q^ or q'^q'( . (For symmetry reasons 
odd powers of qi vanish at midrapidity for central col- 
lisions between equal nuclei.) Such higher order terms 
will affect the 3D fits of the experimentalist, but have no 
effect on equation ((TT)j) . 

We therefore now turn to full three-dimensional Gaus- 
sian fits. We will see that the above analytic expressions 
are easily generalized for this case. 



B. Three-dimensional Gaussian fit algorithm 

Proceeding as in the previous subsection, we start from 
the general three-dimensional Gaussian ansatz jSJ which 
can be written as 

ln(C(g)-l) = InA - {qlRlWsRl+ql R1) ■ (17) 

If the correlation function C[q^^''^ in bin k has error (7^, 
the error on ln(C— 1) is given as in (|ll|l by 



- 1 



(18) 



We minimize 



N 



x^-E 
fc=i 

by setting 



In - 1) - In A + ^ (gf 



dx^ 



0, 



dx^ 



= 



o, s,0- 



91nA 

This leads to a set of 4 coupled linear equations. 



(19) 
(20) 



are directly calculable from the calculated correlator. 
Note that the constant error (7iD,i of the correlator drops 



/3 



(21) 
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FIG. 1: (Color online) One-dimensional slices of the three- 
dimensional correlation function along the "out" , "side" , and 
"long" directions, for pion pairs with K — 0, calculated from 
the blast-wave parameterization. For a given slice, the un- 
plotted g-components equal 1.25 MeV/c (i.e. the center of 
the first bin). The solid (red) curve is the calculated corre- 
lation function from Eq. Q, the dashed (blue) curve shows 
the same slice of the best 3D Gaussian fit l^, with "HBT 
parameters" calculated from the analytic expressions given in 
Sec. lTTTBl 



where a and P take the values 0, o, s, I. The vectors ap- 
pearing here are 



P = (In X,RlRiRf), 
^ _^ln(C(gW)-l) 

k=l 



(22) 
(23) 



V. = +E7T^-l-(C(9^'Vl). (24) 

k=i y'^k) 

while the symmetric 4x4 matrix T has components 

N 



P<Aa> — 



(25) 



In Equations H24|) and 125|) i, j = o,s,l as usual. Note the 
correspondences Vq, ^ X^.i and T^p <-> Y^ i between the 
3D and ID cases. 
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FIG. 2: (Color online) Solid (red) curves show one- 
dimensional slices of the three-dimensional correlation func- 
tion calculated with Eq. Q from the blast-wave parameteri- 
zation, for midrapidity pions with Kt = 0.3 GeV/c. Dashed 
(red) curves show slices of the three-dimensional Gaussian 
form of Equation @, with "HBT parameters" calculated from 
the analytic expressions given in Sec. IIII 51 



The set of linear equations (|21() is easily solved alge- 
braically by diagonalizing the matrix Tap. 



IV. APPLICATION TO BLAST- WAVE MODEL 

Many variants of "hydrodynamically-inspired" models 
of freeze-out have recently been used to calculate spatial 
RMS variances which then were compared to experimen- 
tal HBT radii. A recent example is reported in refer- 
ence The model itself is very simplistic and ignores, 
for example, resonance decay contributions which may be 
important ,2Q] . We ignore such issues with the model it- 
self and simply use it here to discuss differences between 
RMS variances and Gaussian HBT radii. 

We use "realistic" model parameters which best de- 
scribe the data Specifically, we take i?=13.3 fm 
for the source radius, T~97 MeV for the tempera- 
ture, po = l-03 for the maximum transverse flow rapid- 
ity, T = 9 fm/c for the average freeze-out time, and 
Ar = 2.83 fm/c for the emission duration (see [23| for 
details). 



6 




0.05 



0.05 0.05 



"0 0.05 



FIG. 3: (Color online) From the blast-wave parameteriza- 
tion, one-dimensional HBT fit parameters Riu.i and AiD,i 
are calculated with Eqs. 1131141 and plotted as a function of 
the maximum allowed value of any g-component; see text for 
details. Each curve corresponds to one of ten values of Kt- 
0.0, 0.1, 0.2, 0.9 GeV/c. Curves corresponding to high 
Kt are at low (high) values of -RiD.i (Aio.i)- 



S 9 



7 - 



5 
4 

3F 



side 

* Variance 
— HBT radius 

• STAR data 



F 



0.5 0.5 



long 



0.5 




Kj^(GeV/c) 



0.5 

K^(GeV/c) 



FIG. 4: (Color online) One-dimensional HBT fit parameters 
RiD,i and AiD,i as a function of Kt, calculated from the blast- 
wave parametrization with Eqs. 113I14II . For a given Kt, 
the vertical red line represents the variation with fit range 
(see Figure O. Blue stars represent the corresponding radius 
parameters calculated from the RMS variances using Eq. ||SJ . 
Black circles show STAR data with error bars removed 
for clarity. 

A. Correlation functions and analytic fits: results 




side 



0.05 



long 



0.05 




0.05 

q (GeV/c) 




0.05 



FIG. 5: (Color online) From the blast-wave parameterization, 
three-dimensional HBT fit parameters Ri and A are calculated 
with Eqs. 1211 and plotted as a function of the maximum 
allowed value of any g-component; see text for details. Each 
curve corresponds to one of ten values of Kt- 0.0, 0.1, 0.2, 
. . . , 0.9 GeV/c. Curves corresponding to high Kt are at low 
(high) values of Ri (A). The Ri curve for Kt =0 falls above 
the plotting range. 
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FIG. 6: (Color online) Three-dimensional HBT fit parame- 
ters i?iD,i and AiD.i as a function of Kt, calculated from the 
blast- wave parametrization with Eqs. 12111 . For a given Kt, 
the vertical red line represents the variation with fit range 
(see Figure I^J. Blue stars represent the corresponding radius 
parameters calculated from the RMS variances using Eq. JSJ. 
Black circles show STAR data 4], with error bars removed 
for clarity. 



Equation (12) of 22] gives the functional form for the 
single-pion emission function in the blast-wave model. 
Using this for S{K, x), we calculate the correlation func- 
tion for pion pairs with longitudinal pair momentum 
Kl — 0, using a Monte Carlo technique to numerically 
perform the integrals in Eq. (PJ. 

As with experimental data, the correlation function 
is evaluated in finite-sized three-dimensional bins in 
{qo,Qs,qi) of width 2.5 MeV/c in each direction. One- 
dimensional slices of the correlation function in the 
"out" , "side" , and "long" directions are shown in Fig- 



ures ^ and [21 for midrapidity pion pairs with Kt — and 
Kt — 0-3 GeV/ c, respectively. 

The slices of the correlation functions appear quite 
Gaussian, and they are tracked well by the three- 
dimensional Gaussian fit; the fitted correlation strength A 
is very close to 1. The radius parameters calculated from 
the RMS variances © agree quite well with the HBT 
radii extracted from the three-dimensional Gaussian fit 
by solving Eqs. H21|l : both sets are given in the Figures. 
Upon closer inspection one notices, however, that the fit- 
ted outward and longitudinal radii, Ro and especially Ri, 
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tend to be systematically smaller than those extracted 
from the spatial RMS variances; the opposite is true for 
the sideward radii Rs for which the RMS variances give 
slightly smaller values than the Gaussian fit. While these 
differences are small for the blast-wave model parameter- 
ization (at least with the "realistic" parameters studied 
here) , they will be significantly larger (with the same ba- 
sic tendencies as found here) for the hydrodynamic model 
source studied in Sec. Ivl 

The Gaussian fit parameters given in Figs.^andlSlcor- 
respond to using the largest possible q-range in the sums 
over k in Eqs. H24I25|I . discarding only those data points 
for which C is so close to 1 that the Monte Carlo integra- 
tion sometimes yields negative values for C— 1. Due to 
small but noticeable deviations of the correlation func- 
tion from a pure Gaussian, the Gaussian fit parameters 
depend on the number of data points used. We study this 
sensitivity to the fit range in the following subsection. 



B. Fit- range study 

Since no measured correlation function is ever per- 
fectly Gaussian, experimentalists typically perform so- 
called "fit range studies." Here, the measured correlation 
function is fitted with the Gaussian form jSJ, using data 
points in a restricted range of q. With correlation func- 
tions in the one-dimensional quantity Qinv it is common 
to study the variation of fit parameters as the first few 
(lowest-Qinv) data points are left out of the fit. This is 
because statistical fluctuations in these bins may be quite 
large, and due to the visible non-Gaussian nature of the 
measured correlation function there. Three-dimensional 
correlation functions do not suffer from these issues, and 
so usually the experimentalist includes all data points 
with \qi \ < (/max and studies variations of the fit param- 
eters as ^max is varied; any such variations are typically 
folded into systematic errors on the HBT radii. 

Here, we follow the experimentalists' approach. Using 
the correlation function generated from the blast-wave 
model, we calculate HBT parameters from ID and 3D 
Gaussian fits as discussed in Sections IHI Al and IIII Bl re- 
stricting the fc-sums in Eqs. (fT^ . ifT^ . if^ . and if^ 
to include only those data points where all three q- 
components have magnitudes less than Qmax H^- Thus, 
we will not calculate unique HBT radii, but a finite range 
for each fit parameter. 

For various values of Kt, Figure |21 shows the evolution 
of the ID radii with Qmax- Except for Ri at low Kt, the 
parameter variation with fit range is quite mild, corre- 
sponding to a small "non-Gaussian systematic error" on 
the radii. In Figure 01 the range of this variation, indi- 
cated by vertical lines, is plotted as a function of Kt- 
Consistent with the theorem Q that the spatial RMS 
variances ^ of the source control the curvature of the 
correlator C{q) at q = 0, the blue stars in Fig. 0] coincide 
with the (/max limit of the fitted ID radii. The largest 
fit-range variations, indicating the biggest non-Gaussian 



effects in the correlator, are seen at small pair momentum 
Kt- The fit-range sensitivity is most pronounced for Ri 
(where at low Kt it can exceed 0.5 fm) but almost negli- 
gible for Ro and Rg. In short, the ID Gaussian fits to the 
two transverse projections of the correlation function give 
length scales consistent with the spatial RMS variances of 
the source distribution, but non-Gaussian features along 
the longitudinal projection cause the RMS variances to 
overestimate the longitudinal ID HBT radius Ri by up 
to 0.5 fm at low Kt if a reasonable fit range (/max is used 
to extract the latter. This discrepancy is significantly 
larger than the combined statistical and systematic error 
on the experimental value for Ri 4J. 

Figures [S] and |H1 show the same study for the three- 
dimensional Gaussian fits. For reasons explained in 
Sec. IIII AI the non-Gaussian effects in a unified 3D Gaus- 
sian fit are expected to differ from those in ID fits. In- 
deed, in the unified 3D fit non-Gaussian influences also 
appear in Ro, and both Ro and Ri now show fit-range 
variations which exceed the combined statistical and sys- 
tematic errors of the data . The largest fit-range sen- 
sitivity is still seen in the longitudinal direction. In 
Ref. [23| the blast-wave model parameters were deter- 
mined by comparing RMS variances with the measured 
HBT radn (see Figs. ^ and 0), using the experimental 
errors on the latter to extract error estimates for the 
model parameters. The results presented here suggest 
that if the authors had instead compared the measured 
data with HBT radii extracted from a 3D Gaussian fit 
to the calculated correlation function, they would have 
found somewhat different model parameters whose mean 
values in some cases might even have fallen outside the 
likely parameter range quoted in Table II of Ref. . In 
particular, such an "apples-to-apples" comparison may 
allow for somewhat larger fireball lifetimes r and/or emis- 
sion durations At than quoted in Ref. While such 
an improved blast-wave model fit is numerically expen- 
sive and outside the scope of the present paper, it may 
be a worthwhile future project. 



V. HBT RADII FROM HYDRODYNAMICS 

Non- viscous ("ideal") hydrodynamical calculations 
have successfully reproduced differential momentum 
spectra (at least perpendicular to the beam direction) 
at RHIC, including their anisotropics in non-central col- 
lisions and the dependence of these anisotropies on the 
masses of the emitted hadrons As in the blast-wave 
model calculations, very strong collective flow is a criti- 
cal ingredient to reproduce the data. (Of course, in the 
blast-wave parameterization such flow is put in by hand 
while it arises naturally in the hydrodynamical model.) 

Most (but not all [llQ0[l3) hydrodynamic pre- 
dictions of HBT radius parameters have been based on 
calculations of the spatial RMS variances from the hydro- 
dynamically generated emission function, using Eqs. © 
[glllOjI. In spite of the hydrodynamic model's impressive 
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FIG. 7: (Color online) Solid (red) curves show one- 
dimensional slices of the three-dimensional correlation func- 
tion calculated with Eq. from the hydrodynamic model 
with CE Equation of State, for midrapidity pions with Kt ~ 
0.3 GeV/c. Dashed (red) curves show slices of the three- 
dimensional Gaussian form of Equation with "HBT pa- 
rameters" calculated from the analytic expressions given in 

sec. inrsi 



success in describing hadron spectra, these predictions 
of HBT radii were a failure: The calculated longitudinal 
radii Ri were too large (although this problem was less 
severe in Hirano and Tsuda's work ^15.]). while the pre- 
dicted sideward radius Rs was too small, and both Rg 
and Ro showed much less dependence on Kt in theory 
than seen in the data. This, together with similar fail- 
ures by other dynamical models (see P] for a review), has 
become known as the "RHIC HBT Puzzle" . 

Various possibilities to explain and correct this fail- 
ure have been suggested. They include a more realistic 
modeling of the final freeze-out stage |23| , exploration of 
fluctuations in the initial state and ambiguities in the hy- 
drodynamic decoupling criterion , viscous effects due 
to incomplete thermalization (i.e. inapplicability of ideal 
fluid dynamics) , different (more Landau- type) initial 
conditions leading to strong longitudinal hydrodynamic 
acceleration 30] , and the use of more realistic or different 
equations of state (EoS) for the expanding matter [sif . 
None of these suggestions, individually or in combina- 
tion, has been convincingly shown to be able to solve the 
HBT puzzle. Motivated by the blast-wave study in the 
preceding section, we therefore explore here one further 
possibility: that previous comparisons of the data with 
hydrodynamic models might have been misleading since 
the RMS variances from hydro-generated sources differ 



significantly from HBT radii extracted from a Gaussian 
parameterization of the correlation function. Indications 
that this is indeed the case have already emerged from 
the work on ID projections of Hirano and Tsuda and 
Kolb 01, and with our new analytic 3D Gaussian fit al- 
gorithm we can improve on their analysis and study this 
question in more detail. 

For our study of HBT radii from the hydrodynamic 
model we use two different sets of emission functions, 
obtained from running the hydrodynamic code with two 
different equations of state (EoS). Both EoS describe 
the quark-gluon plasma (QGP) as a free gas of mass- 
less particles, but they differ in their treatment of the 
late hadronic stage when the fireball has cooled below 
the critical temperature Tc ~ 165 MeV for hadronization. 
The "CE EoS" [H ^ assumes that the hadron reso- 
nance gas remains not only in thermal, but also in chem- 
ical equilibrium until final kinetic freeze-out. This fails to 
reproduce the observed hadron yields which correspond 
to chemical equilibrium at a tem per ature of about 170 
MeV m. The "NCE EoS" [HM IH takes the im- 
mediate decoupling of hadron abundances at Tc into ac- 
count by introducing non- equilibrium c/iemzca/ potentials 
for each hadron species which ensure that the particle 
yields are held fixed as the temperature and density con- 
tinue to decrease. While the CE EoS was used for the 
hydrodynamic model predictions made for RHIC before 
the accelerator turned on and the hadron abundances 
were measured, the NCE EoS is more realistic and has 
been used in most hydrodynamic studies since 2002. We 
here explore emission functions obtained with either EoS. 

Figures I7I11I present ID projections and ID and 3D 
fit results, analogous to those from the previous section, 
for the emission function from hydrodynamic calculations 
using the CE EoS. Figures fl2llt)l show the same for the 
NCE EoS. Several observations are in order. 

As is apparent from Figures [7| and II 21 the best 3D 
Gaussian fits do not fully reproduce the correlation func- 
tion, even though the correlation function projections 
themselves appear rather Gaussian. Clearly, aspects 
of the correlation function not apparent in the one- 
dimensional projections are partially driving the 3D fit. 
Further, it is interesting to note that, while the pro- 
jections in the "side" direction appear the worst repro- 
duced by the fit, the greatest discrepancy between RMS 
variances and HBT radii are in fact in the "out" and 
"long" directions (c.f. Figures ITTI and [TB|) . Both of these 
points emphasize that the three-dimensional correlator 
can contain important information which docs not ap- 
pear in its one-dimensional projections, and thus in the 
one-dimensional fits. Particularly important in this case 
are strong non-Gaussian features in the longitudinal di- 
rection which cause a significant suppression of the corre- 
lation strength parameter A of the 3D Gaussian fit. This 
in turn creates the appearance of a "bad fit" in the side- 
ward direction even though the ID sideward projection 
looks quite Gaussian itself. 

One draws the same conclusion by examining the fit- 
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FIG. 8: (Color online) From the hydrodynamic model with 
CE EoS, one-dimensional HBT fit parameters RiD,i and AiD,i 
are calculated with Eqs. 1131141 and plotted as a function of 
the maximum allowed value of any g-component; see text for 
details. Each curve corresponds to one of ten values of Kt'- 
0.0, 0.1, 0.2, 0.9 GeV/c. Curves corresponding to high 
Kt are at low (high) values of RiB,i (Aid.O- The Ri curves 
for Kt < 0.2 GeV/c fall above the plotting range. 
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FIG. 10: (Color online) From the hydrodynamic model with 
CE EoS, three-dimensional HBT fit parameters Ri and A are 
calculated with Eqs. 12111 and plotted as a function of the 
maximum allowed value of any g-component; see text for de- 
tails. Each curve corresponds to one of ten values of Kt- 
0.0, 0.1, 0.2, 0.9 GeV/c. Curves corresponding to high 
Kt are at low (high) values of Ri (A). The Ri curves for 
Kt < 0.1 GeV/c fall above the plotting range. 
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FIG. 9: (Color online) One-dimensional HBT fit parame- 
ters -RiD.i and AiD,i as a function of Kt, calculated from the 
hydrodynamic model using CE EoS with Eqs. H1.'-{I14^ . For a 
given Kt, the vertical red line represents the variation with fit 
range (see Figure (HJ. Blue stars represent the corresponding 
radius parameters calculated from the RMS variances using 
Eq. Q. Black circles show STAR data 4], with error bars 
removed for clarity. 



1 - 

6 

5 

41* 

3 - 



out 



side 

Variance 
HBT radius 
STAR data 



0.5 0.5 0.5 

K^(GeV/c} 




0.5 

K^(GeV/c) 



FIG. 11: (Color online) Three-dimensional HBT fit param- 
eters Rio,i and AiD.i as a function of Kt, calculated from 
the hydrodynamic model using CE EoS with Eqs. ll2lj- For a 
given Kt, the vertical red line represents the variation with fit 
range (see Figure EU- Blue stars represent the corresponding 
radius parameters calculated from the RMS variances using 
Eq. @. Black circles show STAR data 4], with error bars 
removed for clarity. 



range systeniatics. As mentioned, non-Gaussian effects 
generate a variation of the HBT parameters with gmax- 
As seen in Figures ISl and [TSl fits in the "out" and "side" 
directions produce ID radii and directional Xi parameters 
which vary very little with g„iax; strong fit-range sensi- 
tivity is only seen in the "long" direction where the ID 
projection deviates most strongly from a Gaussian form. 
In the three-dimensional fits, on the other hand (c.f. Fig- 
ures^] and ^J, the strong non-Gaussian features in the 
qi direction now affect all four fit parameters, generating 
strong fit-range sensitivities also for Ro and A. 



There may (and in general will) be other properties of 
the three-dimensional correlation function to which the 
ID projections and their Gaussian fits are not sensitive 
but which affect the 3D Gaussian fit. The extracted val- 
ues for Ro and thus in general depend significantly on 
the detailed conditions under which the Gaussian fit is 
performed. Hence, a meaningful and accurate compari- 
son between models and experimental data requires that 
the Gaussian fit to the theoretical correlation functions 
is done under similar conditions and constraints (e.g. fit 
range) as the in experiment. 
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FIG. 12: (Color online) Solid (red) curves show one- 
dimensional slices of the three-dimensional correlation func- 
tion calculated with Eq. (O from the hydrodynamic model 
with NCE Equation of State, for midrapidity pious with 
Kt = 0.3 GeV/c. Dashed (red) curves show slices of the 
three-dimensional Gaussian form of Equation with "HBT 
parameters" calculated from the analytic expressions given in 
SecETTH 

VI. DISCUSSION AND CONCLUSIONS 

Let us close with some general observations and sum- 
marize our conclusions. 

Except inasmuch as it couples HBT radii in a 3D fit, 
we have not focused here on the A parameter, since com- 
parison to measurements of A is significantly complicated 
by experimental artifacts . This is also the reason why 
tests of consistency between different experiments gen- 
erally compare HBT radii, not A. In all of the ideal- 
ized calculations presented in this report C(|q|=0) = 2, 
so a purely Gaussian correlation function (generated by 
a purely Gaussian source) would yield A = 1 , with no fit- 
range systematics. Indeed, we find that limq^^^^oA = l 
(see e.g. Figure E)) as expected, but that its value de- 
clines as more bins are included in the fit. In experimen- 
tal data, several factors cause A to fall below its nomi- 
nal value of unity. Our calculations confirm the generally 
held folklore that non-Gaussian effects may be important 
to understanding A. 

Of more fundamental interest are the characteristic 
length scales of the emission region. We have seen 
that RMS variances of model-calculated source func- 
tions, which are frequently compared to experimentally 
extracted HBT radii, may systematically differ from "fit- 
ted" HBT radii which characterize the shape of the cor- 



relation function from the same model. Since the latter 
quantity provides the best "apples-to-apples" comparison 
to published experimental data, this can be an important 
observation. 

Previous attempts 0, 0, |23| to estimate the effect in 
hydrodynamical calculations have focused on numerical 
fits to several one-dimensional projections of the calcu- 
lated correlation function. We here presented an ana- 
lytic method to extract these "ID HBT radii" from the 
projections, and further generalized it to the full three- 
dimensional case. The ID projections represent a set 
of zero measure of the full three-dimensional correlation 
function and, as we have seen, may not be sensitive to im- 
portant three-dimensional information. This information 
influences the unified three-dimensional fit to the corre- 
lation function. Since the unified 3D fit most closely 
mimics the procedure of experimentalists, these effects 
are relevant for comparisons between models and data. 

The magnitude of these effects are model dependent. 
The non-Gaussian nature of emission regions in the blast- 
wave parameterization has been noted before |22j | . It was 
shown here to generate only minor deviations from Gaus- 
sian behaviour in the transverse projections of the cor- 
relation function, but the longitudinal projection shows 
significant non-Gaussian features. In a unified 3D Gaus- 
sian fit, non-Gaussian features were seen to generate fit- 
range sensitivities for all four fit-parameters, leading to 
significant downward shifts of both Ri and i?o, especially 
at low Kt, relative to predictions based on the spatial 
RMS variances of the blast-wave source. 

These tendencies were found to be even more strongly 
exhibited by the HBT radii extracted from hydrodynamic 
model sources. The differences between HBT radii ex- 
tracted from 3D Gaussian fits of the correlator and the 
values l|njl calculated from the spatial RMS variances are 
quite significant and thus relevant in considerations of the 
"RHIC HBT puzzle". In particular, for both equations 
of state considered here, the HBT radii in the "out" and 
"long" directions are significantly lower (and closer to the 
data) than the corresponding RMS variances which have 
been the basis of many "puzzle" discussions (c.f. Fig- 
ures ^] and . As in the blast- wave model, these 3D 
Gaussian fit effects seem to be mostly driven by strong 
non-Gaussian features in the longitudinal projection of 
the correlator. Combining improvements of using the 
NCE EoS and the use of HBT radu instead of RMS vari- 
ances brings the hydrodynamic calculations for the longi- 
tudinal radius Ri into fair agreement with the data over 
the entire measured Kt range. A significant improve- 
ment is also seen in the outward direction, but it is mostly 
concentrated at low Kt, and hence the disagreement be- 
tween the rather steep i^T-dependence of the measured 
Ro radii and the much flatter i^T-dependence of the the- 
oretical results is getting worse. The fitted sideward radii 
Rs show practically no deviation from the corresponding 
RMS variances, and the well-known jl3| problem that 
the hydrodynamically predicted values are significantly 
smaller and show much less i^T-dependence than the 
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FIG. 13: (Color online) From the hydro model with NCE 
EoS, one-dimensional HBT fit parameters Rio.i and AiD,i 
are calculated with Eqs. 1131141 and plotted as a function of 
the maximum allowed value of any g-component; see text for 
details. Each curve corresponds to one of ten values of Kt'- 
0.0, 0.1, 0.2, 0.9 GeV/c. Curves corresponding to high 
Kt are at low (high) values of JJid.z (Aid.O- The Ri curves 
for Kt < 0.1 GeV/c fall above the plotting range. 
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FIG. 15: (Color online) From the hydrodynamic model with 
NCE EoS, three-dimensional HBT fit parameters Ri and A 
are calculated with Eqs. 12111 and plotted as a function of 
the maximum allowed value of any g-component; see text for 
details. Each curve corresponds to one of ten values of Kt'- 
0.0, 0.1, 0.2, . . . , 0.9 GeV/c. Curves corresponding to high 
Kt are at low (high) values of Ri (A). The Ri curves for 
Kt < 0.1 GeV/c fall above the plotting range. 
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FIG. 14: (Color online) One-dimensional HBT fit parameters 
RiD,i and AiD.i as a function of Kt, calculated from the hy- 
drodynamic model using NCE EoS with Eqs. I|13I14|I . For a 
given Kt , the vertical red line represents the variation with fit 
range (see Figure EJ. Blue stars represent the corresponding 
radius parameters calculated from the RMS variances using 
Eq. @. Black circles show STAR data 4], with error bars 
removed for clarity. 
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FIG. 16: (Color online) Three-dimensional HBT fit parame- 
ters RiT>,i and AiD.i as a function of Kt, calculated from the 
hydrodynamic model using NCE EoS with Eqs. I|21|l . For a 
given Kt, the vertical red line represents the variation with fit 
range (see Figure [T^ . Blue stars represent the corresponding 
radius parameters calculated from the RMS variances using 
Eq. Q. Black circles show STAR data 4], with error bars 
removed for clarity. 



data is not alleviated by our improved comparison be- 
tween theory and data. 

While the results presented here can not offer a resolu- 
tion of all aspects of the "RHIC HBT Puzzle" , they refo- 
cus our perception of where the most severe problems are 
located. The strong non-Gaussian effects in qi direction 
and the resulting large downward shift of the fitted lon- 
gitudinal radii (as compared to the corresponding RMS 
variances) largely eliminate the discrepancies between 
hydrodynamically predicted and measured Ri values. A 
number of authors have interpreted the smallness of the 



measured Ri values as evidence for a short fireball lifetime 
T/<10fm/c, inconsistent with the 0(15fm/c) lifetimes 
predicted ^] by the hydrodynamic model. The analy- 
sis presented here resolves this problem. On the other 
hand, even when using the properly extracted Gaussian 
fit values for Rg and Ro and after taking into account the 
resulting decrease of Ro at low Kt, the theoretically pre- 
dicted ratio Ro/Rs is still significantly larger than 1 over 
the entire measured Kt interval, in contradiction to the 
data. Furthermore, the decline of both Ro and Rg with 
increasing pair momentum is still much too weak in the 
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model, in spite of the large transverse flow generated by 
the hydrodynamic expansion. These aspects of the HBT 
Puzzle remain serious and must be addressed by other 
theoretical improvements. 

Finally, one should remember that the raw experimen- 
tal correlation functions hardly ever appear very Gaus- 
sian, due to additional distortions by the final state 
Coulomb interactions between the two charged particles. 
Modern methods of extracting the HBT radii from the 
measured correlator include these Coulomb effects self- 
consistently in the fit function leading to more com- 
plicated (numerical) fit algorithms than the analytical 
one presented in Section ITTll Nonetheless, the measured 
HBT radii extracted from such self-consistent 3D fits 
are affected by non-Gaussian structures in the under- 
lying Bose-Einstein correlations in much the same way 



as discussed here for the simpler case of non-interacting 
particles. Thus, while Coulomb interactions should be 
included in future studies, our analysis should provide 
a good estimate of the direction and magnitude of non- 
Gaussian effects in blast-wave and hydrodynamical mod- 
els, and it points out the importance of such effects in 
the comparison of theory to experiment. 
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